Compare datasets

---

library(raster)

library(plyr)
library(dplyr)
library(reshape)

library(ggplot2)
library(ggpmisc)
library(ggpubr)

library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlBu")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "NA" = "#D5D8DC")

Liens utiles:

• https://www.fao.org/3/x0490E/x0490e0k.htm#annex%203.%20background%20on%20physical%20parameters%20used%20in%20evapotranspiration%20computatio

# era5vA <- read.table("era5_AI.txt") %>%
#   mutate(tas = tas-273.15, pr = mpr * 60 * 60 * 24 * 365, rsds = ssrd / (60*60*24) * 1e-6 * 60*60*24) %>% # conversion from J/m2 to J/sm2 to MJ/m2/day  
#   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","rsds","ea","ET0","AI","cat.AI"))
# 
# 
# write.table(era5vA, "era5vA.txt")

era5vA <- read.table("ERA5/era5_AI.txt") %>% 
  mutate(z = z, tas = t2m - 273.15, sfcWind = wind, model = "historical", source = "ERA5") %>%  #pr, ET0 already in mm/year
  select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)

write.table(era5vA, "era5vA.txt")

wdcvA <- read.table("wdc_AI.txt") %>% mutate(z = z, ea = vapr, ET0 = ET0) %>% # srad in MJ/m2/day, pr and ETO already in mm/year
 select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)


write.table(wdcvA, "wdcvA.txt")


cmip6 <- read.table("cmip6_AI.txt")

cmip6vA <- cmip6 %>% subset(period == "1970_2000") %>%
   group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
   dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(spr, na.rm = T), # pr in mm/year
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>% # ET0 in mm/year
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6vA$cat.AI <- cmip6vA$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6vA$cat.AI[which(cmip6vA$ET0 < 400)] <- "Cold"

write.table(cmip6vA, "cmip6vA.txt")

era5vA <- read.table("era5vA.txt")
wdcvA <- read.table("wdcvA.txt")
cmip6vA <- read.table("cmip6vA.txt")

1 Aridity Index

comp <- merge(select(era5vA, c("lon","lat","source","Continent","z","Rn","ET0","AI","cat.AI")), 
                     select(wdcvA, c("lon","lat","source","Continent","z","Rn","ET0","AI","cat.AI")), by = c("lon","lat")) %>%
        merge(select(cmip6vA, c("lon","lat","source","Continent","z","ET0","AI","cat.AI")), by = c("lon","lat"))

calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

#ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void(base_size = 15)+
    theme(legend.position = "left")

#)

calib$same.cat <- NA
calib$cat.ref <- NA

for(i in 1:nrow(calib)){
  loni = calib$lon[i]
  lati = calib$lat[i]
  catrefi <- subset(calib, lon == loni & lat == lati & source == "CMIP6")$cat.AI %>% as.character()
  cati <- calib$cat.AI[i] %>% as.character()
  calib$same.cat[i] <- ifelse(cati == catrefi, "y",
                              ifelse(cati!= catrefi, "n", 
                                     NA))
  calib$cat.ref[i] <- catrefi
}

write.table(calib, "calib.txt")

# table_calib_c <- calib %>% select(c("Continent","source","cat.AI")) %>%
#   reshape2::melt(id.vars = c("Continent","source")) %>%
#   group_by(Continent,source,value) %>%
#   summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
#   ungroup()
# 
# table_calib_c$value[is.na(table_calib_c$value) == T] <- "NA"
# 
# ggplot(table_calib_c)+
#   geom_col(aes(x = Continent, y = perc, fill = value, position = "fill"))+
#   scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
#   theme_minimal()+
#   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
#   facet_grid(rows = vars(source))+
#   labs(x = "Continent", y = "%", fill = "Aridity\ncategory", y = "Proportion of aridity category")
# 


table_calib <- calib %>% select(c("Continent","Name","Acronym","source","cat.AI")) %>%
  reshape2::melt(id.vars = c("Continent","Name","Acronym","source")) %>%
  group_by(Continent,Name,Acronym,source,value) %>%
  summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
  ungroup() %>% filter(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC"))

table_calib$value[is.na(table_calib$value) == T] <- "NA"
table_calib$Continent[which(table_calib$Continent == "EUROPE-AFRICA")] <- "EUROPE"
table_calib$Continent[which(table_calib$Continent == "CENTRAL-AMERICA")] <- "CENTRAL\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "NORTH-AMERICA")] <- "NORTH\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "SOUTH-AMERICA")] <- "SOUTH\nAMERICA"

table_calib$Continent <- factor(table_calib$Continent, levels = c("AFRICA","EUROPE","ASIA","OCEANIA","NORTH\nAMERICA","CENTRAL\nAMERICA","SOUTH\nAMERICA"))

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = perc, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "%", fill = "Aridity\ncategory")

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = count, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
  scale_y_continuous(position = "right")+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "Count", fill = "Aridity\ncategory")

library(ggrepel)

pie_calib <- table_calib %>% group_by(source, value) %>% summarise(gridcells = sum(count)) %>% mutate(percent = round(gridcells/sum(gridcells)*100, 1), lab.y = (rev(cumsum(rev(gridcells)))) - gridcells*0.5)

ggplot(pie_calib)+geom_bar(aes(x = "", y = gridcells, fill = value), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = percent), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

2 Comparison of datasets

vars <- c("lon","lat","Continent","Name","Acronym","tas","pr","sfcWind","Rn","ea","ET0","AI","cat.AI")
merge.vars <- c("lon","lat","Continent","Name","Acronym")
calib.wide <- select(era5vA, vars) %>%
               merge(select(wdcvA, vars), 
                     by = merge.vars) %>%
          merge(select(cmip6vA, vars), 
                by = merge.vars) %>%
          setNames(c(merge.vars,
                     "tas.era5","pr.era5", "sfcWind.era5","Rn.era5","ea.era5","ET0.era5","AI.era5","cat.AI.era5",
                     "tas.wdc","pr.wdc", "sfcWind.wdc","Rn.wdc","ea.wdc","ET0.wdc","AI.wdc","cat.AI.wdc",
                     "tas.cmip6","pr.cmip6", "sfcWind.cmip6","Rn.cmip6","ea.cmip6","ET0.cmip6","AI.cmip6","cat.AI.cmip6"))

2.1 CMIP6 & ERA5 compared to WDC

2.1.1 Maps

In red: when CMIP6 or ERA5 are warmer than Wordlclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.era5-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: ERA5 - WDC")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - WDC")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or ERA5 are wetter than Worldclim

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or ERA5 compared to Worldclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.era5-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

2.1.2 Scatterplot

Representing for scatterplot only tas, pr and sfcWind (Rn need checking in Wdc, ea needs checking in CMIP6)

cowplot::plot_grid(
  
  ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"))+
  geom_point(aes( y = sfcWind.era5), col = "white")+
  geom_point(aes(y=sfcWind.cmip6), col = "white")+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_void(base_size=15)+
  theme(legend.position = "bottom"),
  
ggplot(calib.wide,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.4, shape = 1)+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.4, shape = 1)+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.4, shape = 1)+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.4, shape = 1)+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn, MJ/M2/m2, WDC, °C", y = "Rn, MJ/M2/m2", col = "Dataset")+
   scale_color_manual(values = c(colorvec[3], colorvec[9]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.4, shape = 1)+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa", col = "Dataset")+
   scale_color_manual(values = c(colorvec[3], colorvec[9]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.4, shape = 1)+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day", col = "Dataset")+
   scale_color_manual(values = c(colorvec[3], colorvec[9]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 2
)

2.1.3 Aridity category difference

calib.wide$same.cat.era5 <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.era5, "y","n"))
calib.wide$change.cat.era5 <- with(calib.wide, paste(cat.AI.era5, cat.AI.cmip6, sep = " to "))
table(calib.wide$change.cat.era5)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                          989                           14 
##                Arid to Humid            Arid to Hyperarid 
##                            5                           14 
##            Arid to Semi-arid                 Cold to Cold 
##                          511                        10377 
##                Cold to Humid         Dry subhumid to Arid 
##                           62                            6 
##         Dry subhumid to Cold Dry subhumid to Dry subhumid 
##                           10                          169 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          382                          146 
##                Humid to Arid                Humid to Cold 
##                            2                          745 
##        Humid to Dry subhumid               Humid to Humid 
##                          219                         5456 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          121                          230 
##       Hyperarid to Hyperarid       Hyperarid to Semi-arid 
##                          725                            3 
##                     NA to NA            Semi-arid to Arid 
##                          570                          114 
##    Semi-arid to Dry subhumid           Semi-arid to Humid 
##                          328                          285 
##       Semi-arid to Semi-arid 
##                          778
calib.wide$dryer.era5 <- with(calib.wide, ifelse(change.cat.era5 == "Arid to Humid", "no",
                                ifelse(change.cat.era5 == "Arid to Semi-arid", "no",
                                ifelse(change.cat.era5 == " Arid to Dry subhumid", "no", 
                                ifelse(change.cat.era5 == "Arid to Hyperarid", "yes",
                                ifelse(change.cat.era5 == "Cold to Humid", "no",
                                ifelse(change.cat.era5 == "Dry subhumid to Cold", "no",
                                ifelse(change.cat.era5 == "Dry subhumid to Humid", "no",
                                ifelse(change.cat.era5 == "Dry subhumid to Arid", "yes",
                                ifelse(change.cat.era5 == "Dry subhumid to Semi-arid", "yes",
                                ifelse(change.cat.era5 == "Humid to arid", "yes",
                                ifelse(change.cat.era5 == "Humid to Semi-arid", "yes",
                                ifelse(change.cat.era5 == "Humid to Cold","no",
                                ifelse(change.cat.era5 == "Humid to Dry subhumid","yes",
                                ifelse(change.cat.era5 == "Semi-arid to Dry subhumid", "no",
                                ifelse(change.cat.era5 == "Semi-arid to Arid", "yes",
                                ifelse(change.cat.era5 == "Semi-arid to Humid", "no",
                                ifelse(change.cat.era5 == "Hyperarid to Arid", "no",
                                ifelse(change.cat.era5 == "Hyperarid to Semi-arid", "no",
                                NA)))))))))))))))))))

ggplot(filter(calib.wide, same.cat.era5 == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.era5))+
  scale_fill_manual(values = c(colorvec[7], colorvec[4]), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 wetter than")+
  ylim(-55,90)+
  theme_void()+
  theme(legend.position = "bottom")

2.2 CMIP6 & WDC compared to ERA5

2.2.1 Maps

In red: when CMIP6 or WDC are warmer than ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.wdc-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: WDC - ERA5")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - ERA5")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or WDC are wetter than ERA5

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.wdc-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or WDC compared to ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.wdc-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

2.2.2 Scatterplot

Representing for scatterplot only tas, pr and sfcWind (Rn need checking in Wdc, ea needs checking in CMIP6)

cowplot::plot_grid(
  
ggplot(calib.wide,aes(x = sfcWind.era5))+geom_point(aes( y = sfcWind.wdc, col = "WDC"))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"))+
  geom_point(aes( y = sfcWind.wdc), col = "white")+
  geom_point(aes(y=sfcWind.cmip6), col = "white")+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_void(base_size=15)+
  theme(legend.position = "bottom"),
  
ggplot(calib.wide,aes(x = tas.era5))+geom_point(aes(y = tas.wdc, col = "WDC"), alpha =0.4, shape = 1)+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
  stat_poly_eq(aes(y = tas.wdc, col = "WDC")) +
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.era5))+geom_point(aes( y = pr.wdc, col = "WDC"), alpha =0.4, shape = 1)+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
    stat_poly_eq(aes(y = pr.wdc, col = "WDC")) +
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.era5))+geom_point(aes( y = sfcWind.wdc, col = "WDC"), alpha =0.4, shape = 1)+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
    stat_poly_eq(aes(y = sfcWind.wdc, col = "WDC")) +
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s", col = "Dataset")+
  scale_color_manual(values = c(colorvec[3], colorvec[9]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.era5))+geom_point(aes( y = Rn.wdc, col = "WDC"), alpha =0.4, shape = 1)+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
    stat_poly_eq(aes(y = Rn.wdc, col = "WDC")) +
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn, MJ/m2/day, ERA5, °C", y = "Rn, W/MJ/m2/day", col = "Dataset")+
   scale_color_manual(values = c(colorvec[3], colorvec[9]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.era5))+geom_point(aes( y = ea.wdc, col = "WDC"), alpha =0.4, shape = 1)+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
    stat_poly_eq(aes(y = ea.wdc, col = "WDC")) +
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa", col = "Dataset")+
   scale_color_manual(values = c(colorvec[3], colorvec[9]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.era5))+geom_point(aes(y = ET0.wdc, col = "WDC"), alpha =0.4, shape = 1)+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha = 0.4, shape = 1)+
    stat_poly_eq(aes(y = ET0.wdc, col = "WDC")) +
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day", col = "Dataset")+
   scale_color_manual(values = c(colorvec[3], colorvec[9]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 2
)

2.2.3 Aridity category difference

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.wdc, cat.AI.cmip6, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid                 Arid to Cold 
##                         1099                           25 
##         Arid to Dry subhumid                Arid to Humid 
##                           36                           77 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           85                          520 
##                 Cold to Cold                Cold to Humid 
##                        10491                          138 
##            Cold to Hyperarid                   Cold to NA 
##                            2                          570 
##         Dry subhumid to Arid         Dry subhumid to Cold 
##                            9                            5 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                          147                          427 
##    Dry subhumid to Semi-arid                Humid to Cold 
##                          146                          588 
##        Humid to Dry subhumid               Humid to Humid 
##                          196                         5112 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          111                           83 
##    Hyperarid to Dry subhumid           Hyperarid to Humid 
##                            1                            1 
##       Hyperarid to Hyperarid       Hyperarid to Semi-arid 
##                          652                            7 
##            Semi-arid to Arid            Semi-arid to Cold 
##                          150                           23 
##    Semi-arid to Dry subhumid           Semi-arid to Humid 
##                          350                          435 
##       Semi-arid to Semi-arid 
##                          775
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "no",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "no",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "no", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "yes",
                                ifelse(change.cat.wdc == "Arid to Cold", "no",       
                                ifelse(change.cat.wdc == "Cold to Humid", "no",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "yes",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "no",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "no",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "yes",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "yes",
                                ifelse(change.cat.wdc == "Humid to arid", "yes",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "yes",
                                ifelse(change.cat.wdc == "Humid to Cold","no",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","yes",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "no",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "yes",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "no",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "no",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "no",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "no",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "no",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "no",
                                NA))))))))))))))))))))))))


ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c(colorvec[7], colorvec[4]), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 wetter than")+
  ylim(-55,90)+
  theme_void()+
  theme(legend.position = "bottom")

3 NDVI

# Download Gimms NDVI
library(gimms)
library(gganimate)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

#ndvi.ref <- downloadGimms(x= 1981,y= 2000,dsn = "NDVI")

list.ndvi <- list.files(path = "NDVI", pattern = "ndvi.*\\.nc4$", full.names = T)[1:2]

ndvi.rast <- rasterizeGimms(list.ndvi, keep = 0) %>% mean(na.rm = T) %>% projectRaster(land_mask) # rasterize Gimms does the quality control and removes NAs
ndvi.df <- ndvi.rast %>% as.data.frame(xy = T) %>% setNames(c("lon","lat","ndvi")) %>% merge(land_mask.df, by = c("lon","lat")) %>% filter(lm != 0)

write.table(ndvi.df, "NDVI/ndvi.df.txt")

ndvi.df <- read.table("NDVI/ndvi.df.txt") 

ggplot(subset(ndvi.df, lm == 1))+geom_tile(aes(x=lon, y = lat, fill = ndvi))+
  scale_fill_continuous_sequential(palette = "Greens 3", rev = T)+
  ylim(-55,90)+labs(fill = "NDVI")+
  theme_void()+theme(legend.position = "bottom")

calibn <- merge(calib, ndvi.df, by = c("lon","lat"))

ggplot(subset(calibn, !is.na(cat.AI)))+geom_boxplot(aes(x=cat.AI, y = ndvi, col = cat.AI))+facet_grid(rows = vars(source))+
  scale_color_manual(values = col.cat)+
  labs(x = "", y = "NDVI")+
  theme_minimal()+theme(legend.position = "none")

aov(ndvi~cat.AI, data = group_by(calibn, source)) %>% summary()
##                Df Sum Sq Mean Sq F value Pr(>F)    
## cat.AI          5   1026  205.17    8708 <2e-16 ***
## Residuals   43200   1018    0.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23577 observations effacées parce que manquantes

range.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(range.ndvi = paste(round(range(ndvi, na.rm = T),2), collapse = " to ")) %>%
  cast(source~cat.AI)

cor.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(pearson.r = round(cor(ndvi, AI, use = "pairwise.complete"),2)) %>%
  cast(source~cat.AI)

knitr::kable(range.ndvi, caption = "NDVI range by source and aridity category")

NDVI range by source and aridity category

source	Cold	Hyperarid	Arid	Semi-arid	Dry subhumid	Humid
CMIP6	-0.13 to 0.77	0.04 to 0.22	0.04 to 0.54	0.01 to 0.82	0.01 to 0.82	-0.01 to 0.92
ERA5	-0.13 to 0.81	0.04 to 0.53	0.01 to 0.46	-0.01 to 0.73	0.04 to 0.82	-0.01 to 0.92
Worldclim	-0.13 to 0.77	0.04 to 0.53	-0.03 to 0.41	-0.04 to 0.72	0.03 to 0.72	-0.08 to 0.92

knitr::kable(cor.ndvi, caption = "Correlation between NDVI and AI by source and aridity category")

Correlation between NDVI and AI by source and aridity category

source	Cold	Hyperarid	Arid	Semi-arid	Dry subhumid	Humid
CMIP6	-0.08	0.15	0.36	0.38	-0.01	0.31
ERA5	-0.16	0.08	0.47	0.30	0.08	0.38
Worldclim	0.05	0.05	0.49	0.39	0.14	0.52